%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Decomposition.m   (Script to generate Table 9)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%--------------------------------------------------------------------------
% 1) Declare parameters
%--------------------------------------------------------------------------

global alphah alphal mh ml Q nu
global zeta lambda eta delta beta delta_p
global zbar chi


%--------------------------------------------------------------------------
% 2) Compute the baseline moments (M0) for 1981–85
%--------------------------------------------------------------------------

M0 = struct();
M0.scope_avg = omega_avg_initial;     % omega_avg for 1981–85
M0.RD_H    = rd_Y_h_initial;           % rd_Y_h for 1981–85
M0.RD_L    = rd_Y_l_initial;           % rd_Y_l for 1981–85
M0.trans   = trans_share_initial;      % trans_share for 1981–85
M0.growth  = g_avg_initial^(zeta/(zeta+lambda)) - 1;  % annual growth for 1981–85


%--------------------------------------------------------------------------
% 3) Loop over key parameters and generate moments
%--------------------------------------------------------------------------

deviations = zeros(6, 5);

for i = 1 : 6
    if i==1
        phi=phi_ending; theta=theta_ending; ss=ss_initial; mu=mu_initial; iota=iota_initial; gamma=gamma_initial; rho=rho_initial;
    elseif i==2
        phi=phi_ending; theta=theta_initial; ss=ss_initial; mu=mu_initial; iota=iota_initial; gamma=gamma_initial; rho=rho_initial;
    elseif i==3
        phi=phi_initial; theta=theta_ending; ss=ss_initial; mu=mu_initial; iota=iota_initial; gamma=gamma_initial; rho=rho_initial;
    elseif i==4
        phi=phi_initial; theta=theta_initial; ss=ss_ending; mu=mu_initial; iota=iota_initial; gamma=gamma_initial; rho=rho_initial;
    elseif i==5
        phi=phi_initial; theta=theta_initial; ss=ss_initial; mu=mu_ending; iota=iota_ending; gamma=gamma_initial; rho=rho_initial;
    else
        phi=phi_initial; theta=theta_initial; ss=ss_initial; mu=mu_initial; iota=iota_initial; gamma=gamma_ending; rho=rho_ending;        
    end

    %--------------------------------------------------------------------------
    %  A) SOLVE THE GE EQUILIBRIUM FOR THOSE PARAMETERS
    %--------------------------------------------------------------------------

    geq_vec=@(input)geq_solver_policy(input,iota,mu,gamma,rho,phi,theta);
    [answer,value,exitflag] = fsolve(geq_vec, guess, options);
    if exitflag ~= 1
        disp('First-Order Optimality is not Achieved')
    end
    
    % Construct some output from the g.e. problem
    ih=answer(1);
    il=answer(2);
    omegah=answer(3);
    omegal=answer(4);
    a=answer(5);
    v1h=answer(6);
    v1l=answer(7);
    ns=answer(8);

    %--------------------------------------------------------------------------
    %  B) RECOMPUTE ALL “INTERMEDIATE” QUANTITIES EXACTLY AS IN MAIN PROGRAM
    %--------------------------------------------------------------------------

    omega_avg=alphah*omegah+alphal*omegal;
    nbh=alphah*(1-ih*X(omegah));
    nbl=alphal*(1-il*X(omegal));
    nb=nbh+nbl;
    sigmah=nbh/nb;
    sigmal=nbl/nb;
    mb=phi*(ns/nb)^nu;
    ms=phi*(nb/ns)^(1-nu);
    
    eh=(alphah*mh*a+alphal*ml)*(alphah*Q(1,1)+alphal*Q(2,1))/((alphah*mh+alphal*ml)*(alphah*Q(1,1)*a+alphal*Q(2,1)));
    el=(alphah*mh*a+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2))/((alphah*mh+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2)*a));
    gh=1+gamma*eh*(ih*X(omegah)+(1-ih*X(omegah))*mb);
    gl=1+gamma*el*(il*X(omegal)+(1-il*X(omegal))*mb);
    
    g_avg=gh*(alphah*Q(1,1)+alphal*Q(2,1)/a)/(alphah*Q(1,1)+alphal*Q(2,1));
    rm=beta/(g_avg^(epsilon*zeta/(lambda+zeta)));
    rt=1/rm-1+delta;
    w=lambda*((eta/rt)^(eta/(zeta+lambda)))*((alphah*mh+alphal*ml)*g_avg*zbar)^(zeta/(lambda+zeta));
    
    o=(ms*theta*(sigmah*v1h+sigmal*v1l)*gamma)/(1-(1-ms*theta)*rm*delta_p*g_avg^(zeta/(lambda+zeta)));
    
    zlbar=(alphah*mh+alphal*ml)*zbar/(alphah*mh*a+alphal*ml);
    zhbar=a*zlbar;
        
    rd_exp_h=(1-ss)*chi*(zbar^(zeta/(zeta+lambda)))*(ih^(1+rho))/(1+rho);
    rd_exp_l=(1-ss)*(zbar^(zeta/(zeta+lambda)))*(il^(1+rho))/(1+rho);
    
    zhr=(alphah*Q(1,1)*zhbar+alphal*Q(2,1)*zlbar)/(alphah*Q(1,1)+alphal*Q(2,1));
    zlr=(alphah*Q(1,2)*zhbar+alphal*Q(2,2)*zlbar)/(alphah*Q(1,2)+alphal*Q(2,2));
    
    Yh=mh*gh*zhr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));
    Yl=ml*gl*zlr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));
    
    rd_Y_h=rd_exp_h/Yh;
    rd_Y_l=rd_exp_l/Yl;
 
    patent_count=alphah*ih+alphal*il;
    trans_share=(alphah*ih*(1-X(omegah))+alphal*il*(1-X(omegal)))*(ms*(1-((1-ms)*delta_p)^10))/((1-(1-ms)*delta_p)*patent_count);

    if v1h*gamma >= o && v1l*gamma>=o
        disp('Both types keep patents');
    elseif v1h*gamma > o && v1l*gamma<o
        disp('High type keeps patents; low type sells patents');
    else
        disp('Both types sell patents');
    end

    % Package the five moments:
    M1 = struct();
    M1.scope_avg = omega_avg;
    M1.RD_H    = rd_Y_h;
    M1.RD_L    = rd_Y_l;
    M1.trans   = trans_share;
    M1.growth  = g_avg^(zeta/(zeta+lambda))-1;

    % Compute percentage deviations (or percentage‐point for growth):
    deviations(i,1) = (M1.scope_avg - M0.scope_avg) / (range_avg2-range_avg) * 100;   % %Δ scope_H
    deviations(i,2) = (M1.RD_H    - M0.RD_H   ) / (rd_sales_h2-rd_sales_h)    * 100;   % %Δ RD_H
    deviations(i,3) = (M1.RD_L    - M0.RD_L   ) / (rd_sales_l2-rd_sales_l)     * 100;   % %Δ RD_L
    deviations(i,4) = (M1.trans   - M0.trans  ) / (trans_target2-trans_target)   * 100;   % %Δ trans‐share
    deviations(i,5) = (M1.growth  - M0.growth ) * 100;                % Δ growth (pp)
end


%--------------------------------------------------------------------------
% 4) Generate Table 9
%--------------------------------------------------------------------------

rowNames = {
    'Patent Market (ϕ,θ)';
    'Matching Efficiency (ϕ)';
    'Sellers Bargaining Power (θ)';
    'R&D Tax Credit (σ)';
    'Production Cost (µ,ι)';
    'Rare Good Ideas (γ,ρ)'
};

varNames = {
    'Prod. Scope';
    'R&D(H)';
    'R&D(L)';
    'Patent Trade';
    'Growth'
};

Table9 = array2table(deviations, ...
    'RowNames',    rowNames, ...
    'VariableNames', varNames);

disp(Table9);
writetable(Table9,  "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table9.xlsx", 'WriteRowNames', true);



